#Calculate across-depths average regional R1
R1 <- do.call(rbind, myelin.glasser.7T) #R1 at analyzed depths, 20%-80% of thickness
R1 <- R1 %>% select(contains("ROI"))
R1 <- colMeans(R1) %>% as.data.frame()  %>% set_names("mean.R1")
R1 <- R1 %>% mutate(orig_parcelname = row.names(R1))

#Save out lh and rh mean R1 data 
R1.lh <- R1$mean.R1[1:180]
write.table(R1.lh, "/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/R1.lh.glasser.txt", col.names = F, row.names = F, quote = F)
R1.rh <- R1$mean.R1[181:360]
write.table(R1.rh, "/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/R1.rh.glasser.txt", col.names = F, row.names = F, quote = F)

df.list <- list(R1, T1wT2w, MTsat, MBP) #dfs to merge
myelin.maps <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
myelin.maps <- merge(myelin.maps, glasser.regions, by = c("orig_parcelname"), sort = F)
myelin.maps <- myelin.maps[!(myelin.maps$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/R1_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.9, height = 2.7)

#run PCA
maps.PCA <- myelin.maps %>% select(label, T1wT2w, MTsat, MBP)
myelin.map.PCA <- maps.PCA %>% select(-label) %>% prcomp(., scale = T, center = T)

#variance explained
summary(myelin.map.PCA)

#PC1
myelin.map.PC1 <- myelin.map.PCA$x[,1] %>% as.data.frame %>% purrr::set_names("PC1")
myelin.map.PC1$label <- myelin.maps$label
myelin.map.PC1$PC1 <- myelin.map.PC1$PC1*-1 #sign of loadings is arbitrary; flip for positive correlation

myelin.map.PC1 %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = PC1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-1.2, 3.5), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/PC1_myelinmaps_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.9, height = 2.7)

myelin.maps <- merge(myelin.maps, myelin.map.PC1, sort = F)

write.csv(myelin.maps, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/myelin_maps.csv", quote = F, row.names = F)

cor.test(myelin.maps$mean.R1, myelin.maps$PC1, method = c("spearman"))

myelin.maps %>%
ggplot(., aes(x = mean.R1, y = PC1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nR1", y="Myelin Map PC1\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/PC1_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.65, height = 1.6)

myelin.maps.regression <- lm(scale(myelin.maps$PC1) ~ scale(myelin.maps$mean.R1))
summary(myelin.maps.regression)$r.squared

#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
  meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
  meanmap$orig_parcelname <- row.names(meanmap)
  meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
  return(meanmap)
}

regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
  regional_means(depth)
})

regional.myelin.alldepths <-  do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- as.numeric(substr(row.names(regional.myelin.alldepths), 1, 1)) #assign depth variable

regional.myelin.alldepths <-  regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

area4.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "4_ROI")
area46.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "46_ROI")
areaa24.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "a24_ROI")

examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)

ggplot() +
  geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = depth, group = region, color = region), se = F, linewidth = .8) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  theme_classic() +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Area_profiles_group.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.55)

#R1 values in all regions at all depths for all participants
myelin.glasser.7T.alldepths <-  do.call(rbind, myelin.glasser.7T)
myelin.glasser.7T.alldepths$depth <- as.numeric(substr(row.names(myelin.glasser.7T.alldepths), 1, 1)) #assign depth variable

#Wide to long format
cols_to_pivot <- names(myelin.glasser.7T.alldepths)[grepl("ROI", names(myelin.glasser.7T.alldepths))]
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1")

#Select exemplar regions of interest
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI" | orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI" | orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")

#Format grouping variables
myelin.glasser.7T.alldepths$subject_id <- as.factor(myelin.glasser.7T.alldepths$subject_id)
myelin.glasser.7T.alldepths$orig_parcelname <- as.factor(myelin.glasser.7T.alldepths$orig_parcelname)

#Plot for the oldest 3 participants 
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(age > 32)

ggplot() +
  geom_smooth(data = myelin.glasser.7T.alldepths, aes(x = R1, y = depth, color = orig_parcelname, group = interaction(orig_parcelname, subject_id)), se = F, linewidth = .4) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Area_profiles_individual.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

#R1 values in all regions at all depths for all participants
myelin.glasser.micapni.alldepths <- do.call(rbind, myelin.glasser.micapni)
regional.myelin.alldepths.micapni <- myelin.glasser.micapni.alldepths %>% group_by(depth) %>% summarise(across(1:360, \(x) mean(x, na.rm = TRUE))) #average depth R1 per region 
regional.myelin.alldepths.micapni <- pivot_longer(regional.myelin.alldepths.micapni, cols = 2:361, names_to = "label", values_to = "mean.R1")
regional.myelin.alldepths.micapni <- left_join(regional.myelin.alldepths.micapni, glasser.regions, by = "label")

regional.myelin.alldepths.micapni <-  regional.myelin.alldepths.micapni[!(regional.myelin.alldepths.micapni$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

area4.myelin.alldepths <- regional.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "4_ROI")
area46.myelin.alldepths <- regional.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "46_ROI")
areaa24.myelin.alldepths <- regional.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "a24_ROI")

examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)

ggplot() +
  geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = as.numeric(depth), group = region, color= region), se = F, linewidth = .8) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  scale_x_continuous(breaks = c(0.5, 0.55, 0.6, 0.65)) +
  scale_y_reverse(breaks = c(2, 4, 6, 8, 10, 12,  14), limits = c(14, 2)) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Area_profiles_group_micapni.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.55)

R1.ribbon.cov <- regional.myelin.alldepths %>% group_by(label) %>% do(profile.cov.7T = EnvStats::cv(.$mean.R1)) %>% unnest(cols = c("profile.cov.7T"))

R1.ribbon.cov.micapni <- regional.myelin.alldepths.micapni %>% group_by(label) %>% do(profile.cov.micapni = EnvStats::cv(.$mean.R1)) %>% unnest(cols = c("profile.cov.micapni"))

R1.ribbon.cov <- left_join(R1.ribbon.cov, R1.ribbon.cov.micapni, by = "label")
write.csv(R1.ribbon.cov, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/R1_cov.csv", quote = F, row.names = F)

cor.test(R1.ribbon.cov$profile.cov.7T, R1.ribbon.cov$profile.cov.micapni, method = "spearman")

R1.ribbon.skewness <- regional.myelin.alldepths %>% group_by(label) %>% do(profile.skewness.7T = skewness(.$mean.R1)) %>% unnest(cols = c("profile.skewness.7T"))

R1.ribbon.skewness.micapni <- regional.myelin.alldepths.micapni %>% group_by(label) %>% do(profile.skewness.micapni = skewness(.$mean.R1)) %>% unnest(cols = c("profile.skewness.micapni"))

R1.ribbon.skewness <- left_join(R1.ribbon.skewness, R1.ribbon.skewness.micapni, by = "label")
write.csv(R1.ribbon.skewness, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/R1_skewness.csv", quote = F, row.names = F)

cor.test(R1.ribbon.skewness$profile.skewness.7T, R1.ribbon.skewness$profile.skewness.micapni, method = "spearman")

R1.ribbon.skewness %>%
ggplot(., aes(x = profile.skewness.7T, y = profile.skewness.micapni)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.5, 1.5), breaks = c(-0.5, 0, 0.5, 1, 1.5)) +
  scale_y_continuous(limits = c(-0.5, 1.5), breaks = c(-0.5, 0, 0.5, 1, 1.5)) +
  labs(x="\nSkewness: Developmental", y="Skewness: High-resolution\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Skewness_correlationplot.pdf", device = "pdf", width = 2.68, height = 1.6)

source("/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/surface_metrics/surface_measures/extract_surfacestats.R") 
CT.glasser <- extract_surfacestats(myatlas = "glasser", mymeasure = "ThickAvg") #cortical thickness in glasser regions for each scan in the final study sample, generated by fs-tabulate pipeline
CT.glasser.frontal <- CT.glasser[, names(CT.glasser) %in% glasser.frontal$orig_parcelname]
CT.glasser.frontal <- CT.glasser.frontal[, !(names(CT.glasser.frontal) %in% glasser.snr.exclude$orig_parcelname)]
CT.glasser.frontal$subject_id <- CT.glasser$subject_id
CT.glasser.frontal$session_id <- CT.glasser$session_id
CT.glasser.frontal$age <- CT.glasser$age

average.CT.frontal <- CT.glasser.frontal %>% select(contains("ROI")) %>% colMeans() %>% cbind() %>% as.data.frame() %>% purrr::set_names("CT") %>% mutate(orig_parcelname = row.names(.)) %>% merge(., glasser.regions)

mean(average.CT.frontal$CT)

sd(average.CT.frontal$CT)

ggplot(average.CT.frontal, aes(x = CT)) +
  geom_density(linewidth = 0.5) +
  theme_classic() +
  labs(x="\nCortical Thickness", y="Density\n") +
  theme(
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

region_cols <- setdiff(names(CT.glasser.frontal), c("subject_id", "session_id", "age"))
age.breaks <- seq(10, 32.5, length.out = 6)

average.CT.frontal.agebinned <- CT.glasser.frontal %>%
  mutate(age_bin = cut(age, breaks = age.breaks, include.lowest = TRUE)) %>% 
  group_by(age_bin) %>%
  summarise(
    mean.age = mean(age, na.rm = TRUE),
    n = dplyr::n(),
    across(all_of(region_cols), ~ mean(.x, na.rm = TRUE)),
    .groups = "drop"
)

average.CT.frontal.agebinned <- average.CT.frontal.agebinned %>% pivot_longer(-c(age_bin, mean.age, n), names_to = "region", values_to = "CT")

bin1 <- average.CT.frontal.agebinned %>% filter(age_bin == "[10,14.5]")
bin2 <- average.CT.frontal.agebinned %>% filter(age_bin == "(14.5,19]")
bin3 <- average.CT.frontal.agebinned %>% filter(age_bin == "(19,23.5]")
bin4 <- average.CT.frontal.agebinned %>% filter(age_bin == "(23.5,28]")
bin5 <- average.CT.frontal.agebinned %>% filter(age_bin == "(28,32.5]")

ggplot() +
  geom_density(data = bin5, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[5] ) +
  geom_density(data = bin4, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[4]) +
  geom_density(data = bin3, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[3]) +
  geom_density(data = bin2, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[2]) +
  geom_density(data = bin1, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[1]) +
  theme_classic() +
  labs(x="\nCortical thickness (mm)", y="Density\n") +
  theme(
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Frontal_CT_agebins.pdf", device = "pdf", width = 2.1, height = 1.8)

ribbon_voxelcounts <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/ribbon_voxelcounts/ribbon_FS_voxelcounts.RDS") #voxel count distributions in the frontal lobe for each longitudinal study session
ribbon_voxelcounts <- left_join(participants, ribbon_voxelcounts, by = c("subject_id", "session_id")) #now for only our final study sample

#Identify the most frequent number of compartments in the cortical ribbon for all study participants
ribbon.compartmentnum.mode <- ribbon_voxelcounts %>% group_by(subject_id, session_id) %>% do(compartment.number = .$voxelcount[which.max(.$count.percents)]) %>% unnest(cols = c("compartment.number"))

table(ribbon.compartmentnum.mode$compartment.number)

#Calculate the percent of vertices that show each number of unique voxel compartments
average.counts.frontal <-  ribbon_voxelcounts %>% group_by(voxelcount) %>% do(mean.percent = mean(.$count.percents)) %>% unnest(cols = c("mean.percent"))

kable(average.counts.frontal)

#98% of vertices have 2 or more unique voxels
average.counts.frontal %>% filter(voxelcount > 1) %>% colSums(.)

#73% of vertices have 3 or more unique voxels
average.counts.frontal %>% filter(voxelcount > 2) %>% colSums(.)

#Cumulative percentages plot
average.counts.frontal <- average.counts.frontal %>%
  arrange(desc(voxelcount)) %>%
  mutate(cumulative_percentage = cumsum(mean.percent)) %>%
  arrange(voxelcount)

average.counts.frontal %>% filter(voxelcount < 6) %>% #exclude .1% of vertices with 6+ 
ggplot(., aes(x = cumulative_percentage, y = voxelcount, xend = 0)) +
  geom_segment(color = myelin.colorbar[4]) +
  geom_point(size = 1) +
  theme_classic() +
  labs(x="\nCumulative percentage", y="Unique compartments\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Cumulativepercentage_compartments.pdf", device = "pdf", width = 1.7, height = 1.8)

average.counts.frontal %>% filter(voxelcount < 6) %>% #exclude .1% of vertices with 6+ 
ggplot(., aes(x = voxelcount, y = mean.percent)) +
  geom_col(fill = myelin.colorbar[4]) +
  theme_classic() +
  labs(x="\nUnique Compartments", y="Percent of Frontal Vertices\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ribbon_voxelcounts.agebinned <- ribbon_voxelcounts %>% filter(voxelcount < 6) %>% mutate(age_bin = cut(age, breaks = age.breaks, include.lowest = TRUE)) %>% 
  group_by(age_bin, voxelcount) %>% 
  do(mean.age = mean(.$age), mean.percent = mean(.$count.percents)) %>% unnest(cols = c("mean.age", "mean.percent"))

bin1 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "[10,14.5]")
bin2 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(14.5,19]")
bin3 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(19,23.5]")
bin4 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(23.5,28]")
bin5 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(28,32.5]")

ggplot() +
  geom_col(data = bin5, aes(x = voxelcount, y = mean.percent), width = .8, fill = myelin.colorbar[5], se = F ) +
  geom_col(data = bin4, aes(x = voxelcount, y = mean.percent), width = 0.65, fill = alpha(myelin.colorbar[4], 1), se = F) +
  geom_col(data = bin3, aes(x = voxelcount, y = mean.percent), width = 0.5, fill = alpha(myelin.colorbar[3], 1), se = F) +
  geom_col(data = bin2, aes(x = voxelcount, y = mean.percent), width = 0.35, fill = alpha(myelin.colorbar[2], 1), se = F) +
  geom_col(data = bin1, aes(x = voxelcount, y = mean.percent), width = 0.2, fill = alpha(myelin.colorbar[1], 1), se = F) +
  theme_classic() +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5)) +
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50)) +
  labs(x="\nUnique compartments", y="Percent of frontal vertices\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Frontal_compartments_agebins.pdf", device = "pdf", width = 2.1, height = 1.8)

#Average depth-dependent R1 profile for the frontal lobe in high resolution (.5 mm) data
frontal.myelin.alldepths.micapni <- regional.myelin.alldepths.micapni[regional.myelin.alldepths.micapni$orig_parcelname %in% glasser.frontal$orig_parcelname,]
depthR1.micapni <- frontal.myelin.alldepths.micapni %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1"))

ggplot(data = depthR1.micapni, aes(x = depth.R1, y = depth)) +
  geom_point() +
  geom_smooth(linewidth = .8, se = F, method = "gam", formula =y ~ s(x, k = 7), color = myelin.colorbar[4]) +
  theme_classic() +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)) +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

#Fit a gam to the depth-by-R1 data and calculate the second derivative (rate of change in the rate of change to identify transition zones)
depthR1.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.micapni) 
depthR1.derivatives.micapni <- derivatives(depthR1.model.micapni, order = 2, eps = .1, interval = "simultaneous") #second derivative

ggplot() +
  geom_path(data = depthR1.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.3) +
  theme_classic() +
  scale_x_continuous(breaks = c(0, 0.003, 0.006)) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)) +
  xlab("\nSecond derivative") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/R1_secondderivative_micapni.pdf", device = "pdf", width = 1.6, height = 1.8)

#Repeat the above for our exemplar regions to see if the same transition zones emerge
depthR1.4.micapni <- frontal.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI") %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1")) %>% mutate(region = "4_ROI")
depthR1.4.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.4.micapni) 
depthR1.4.derivatives.micapni <- derivatives(depthR1.4.model.micapni, order = 2, eps = .1, interval = "simultaneous") 

depthR1.46.micapni <- frontal.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI") %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1")) %>% mutate(region = "46_ROI")
depthR1.46.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.46.micapni) 
depthR1.46.derivatives.micapni <- derivatives(depthR1.46.model.micapni, order = 2, eps = .1, interval = "simultaneous") 

depthR1.a24.micapni <- frontal.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI") %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1")) %>% mutate(region = "a24_ROI")
depthR1.a24.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.a24.micapni) 
depthR1.a24.derivatives.micapni <- derivatives(depthR1.a24.model.micapni, order = 2, eps = .1, interval = "simultaneous") 

ggplot() +
  geom_path(data = depthR1.a24.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.2, color = myelin.colorbar[2]) +
  geom_path(data = depthR1.46.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.2, color = myelin.colorbar[3]) +
  geom_path(data = depthR1.4.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.2, color = myelin.colorbar[5]) +
  theme_classic() +
  scale_x_continuous(breaks = c(0, 0.003, 0.006)) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)) +
  xlab("\nSecond derivative") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2_supplementary/R1_secondderivative_micapni_exemplarregions.pdf", device = "pdf", width = 2.2, height = 2.4)

#Identify change points in the relationship between 2 variables. Segmented creates broken-line regression models where relationship between x and y is piece-wise linear and broken up via break points. It identifies where there is a break in the relationship between x and y 
depthR1.lm.micapni = lm(depth.R1 ~ depth, data = depthR1.micapni)
break.points <- segmented::segmented(depthR1.lm.micapni, seg.Z = ~depth, npsi = 2) #identify to depth break points via change point analysis
break.points$indexU

#Average depth-dependent R1 profile for the frontal lobe in 1 mm data
frontal.myelin.alldepths <- regional.myelin.alldepths[regional.myelin.alldepths$orig_parcelname %in% glasser.frontal$orig_parcelname,]
depthR1.7T <- frontal.myelin.alldepths %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1"))

ggplot(data = depthR1.7T, aes(x = depth.R1, y = depth)) +
  geom_point() +
  geom_smooth(linewidth = .8, method = "gam", formula = y ~ s(x, k = 6), color = myelin.colorbar[4]) + 
  theme_classic() +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

#Fit a gam to the depth-by-R1 data and calculate the second derivative
depthR1.model <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.7T)
depthR1.derivatives <- derivatives(depthR1.model, order = 2, eps = .1, interval = "simultaneous") #second derivative

ggplot(depthR1.derivatives, aes(x = abs(derivative), y = data)) +
  geom_path(linewidth = 0.3) + 
  theme_classic() +
  scale_x_continuous(breaks = c(0, 0.005, 0.01)) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  xlab("\nSecond derivative") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2_supplementary/R1_secondderivative_7T.pdf", device = "pdf", width = 2.2, height = 2.4)

myelin.smd.7T <- do.call(rbind, myelin.compartments.7T) #R1 in superficial, middle, deep for all participants and regions

myelin.smd.7T <- myelin.smd.7T[ , 
  names(myelin.smd.7T) %in% glasser.frontal$orig_parcelname] #just frontal
myelin.smd.7T <- myelin.smd.7T[ , 
  !(names(myelin.smd.7T) %in% glasser.snr.exclude$orig_parcelname)] #good SNR

myelin.smd.7T$depth <- str_remove(row.names(myelin.smd.7T), "\\..*")

cols_to_pivot <- names(myelin.smd.7T)[grepl("ROI", names(myelin.smd.7T))]
myelin.smd.7T <- myelin.smd.7T %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1") #long format
myelin.smd.7T <- left_join(myelin.smd.7T, glasser.regions, by = "orig_parcelname")

myelin.smd.7T <- myelin.smd.7T %>% group_by(depth, label) %>% do(mean.depth.R1 = mean(.$R1)) %>% unnest(cols = c("mean.depth.R1")) #R1 in superficial, middle, deep in frontal regions averaged across participants

myelin.smd.7T %>% 
  ggplot(aes(x = mean.depth.R1, y = depth)) +
  geom_violin(aes(fill = depth), color = "white", bw = .006) +
  stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
  theme_classic() +
  scale_fill_manual(values = c("#004A38FF", "#809dc4", "#A29DC4")) +
  scale_x_continuous(limits = c(0.497, 0.65), expand = c(0,0)) +
  scale_y_discrete(expand = c(0, 0)) +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Depthplot_regionalR1_frontal.pdf", device = "pdf", dpi = 500, width = 2.2, height = 1.7)

for(this.depth in c(unique(myelin.smd.7T$depth))){
  
  depth.data <- myelin.smd.7T %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
  
  depth.plot <- ggplot() + 
  geom_brain(data = depth.data, atlas = glasser,  mapping = aes(fill = mean.depth.R1), colour=I("gray10"), size=I(.05)) + theme_classic() + theme(legend.position = "none") + 
   paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.48, 0.6), oob = squish, na.value = "white") + theme(legend.position = "none") +
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray10", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 

  print(depth.plot)
  ggsave(filename = sprintf("/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/%s_corticalmap_frontal.pdf", this.depth), device = "pdf", dpi = 500, width = 4.1, height = 1.5)
}